Incremento del nivel medio del mar

Juan Sebastián Mendoza Páez

Estudiante Ing. Forestal

¿Por qué es importante?

NASA

“Desde 1880, el nivel del mar global ha aumentado 20 centimetros. Para el 2100 se proyecta que aumente entre 30 y 122 centimetros más.”

Bases de datos

Fuentes

  • Oficina Nacional de Administración Oceánica y Atmosférica (NOAA)
  • Organización de Investigación Científica e Industrial del Commonwealth (CSIRO)
Mostrar Código
pacman::p_load(TSstudio, tsbox, tsoutliers, tidyverse, equatiomatic,
               lubridate, xts, magrittr, forecast, dygraphs)

url <- 'https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt'
cov <- read.table(file = url) %>%
  select(3, 4) %>%
  rename(fecha = 1, co2 = 2) %>%
  mutate(fecha = date_decimal(fecha) %>%
           floor_date(unit = 'month')) %>%
  filter(fecha <= '2013-12-01')

data <- read.csv('gmsl.csv') %>%
  mutate(fecha = date_decimal(Time) %>%
           floor_date(unit = 'month')) %>%
  rename(gsml = 2) %>%
  select(4, 2)%>%
  filter(fecha > '1958-02-01') %>%
  left_join(y = cov) %>%
  {xts(x = .[,2:3], order.by = .$fecha)}

rm(cov)

dygraph(data, main = 'Nivel medio del mar global vs
        Concentración de CO2 en la atmósfera', height = 400) %>%
  dySeries("co2", axis = 'y2', label = 'CO2') %>%
  dySeries('gsml', label = 'GSML') %>%
  dyAxis('y', label = 'Nivel del mar (mm)') %>%
  dyAxis('y2', label = 'CO2 (ppm)', independentTicks = T) %>%
  dyRangeSelector()

Estacionalidad y prueba de raíces unitarias

Mostrar Código
acf(train[,'gsml'], main = '')

Mostrar Código
tseries::adf.test(x = train[,'gsml'])

    Augmented Dickey-Fuller Test

data:  train[, "gsml"]
Dickey-Fuller = -2.7778, Lag order = 8, p-value = 0.2491
alternative hypothesis: stationary

Modelos

Mostrar Código
matriz_diseño <- model.matrix(~ -1 + train[,'co2'])
colnames(matriz_diseño) <- 'co2'

modelo <- Arima(y = train[,'gsml'], 
                order = c(2,1,3),
                xreg = matriz_diseño,
                include.drift = T)

modelo %>% lmtest::coeftest()

z test of coefficients:

       Estimate Std. Error  z value  Pr(>|z|)    
ar1    0.797070   0.054991  14.4947 < 2.2e-16 ***
ar2   -0.095533   0.054919  -1.7395   0.08194 .  
ma1   -0.339090   0.040237  -8.4274 < 2.2e-16 ***
ma2    0.208116   0.038717   5.3753 7.646e-08 ***
ma3   -0.693544   0.028246 -24.5535 < 2.2e-16 ***
drift  0.208518   0.041471   5.0281 4.955e-07 ***
co2   -0.171228   0.087716  -1.9521   0.05093 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC: 2595.726
Mostrar Código
mod <- Arima(y = train[,'gsml'], 
                order = c(2,1,3),
                include.drift = T)

mod %>% lmtest::coeftest()

z test of coefficients:

       Estimate Std. Error  z value  Pr(>|z|)    
ar1    0.794860   0.055187  14.4031 < 2.2e-16 ***
ar2   -0.092967   0.055205  -1.6840   0.09218 .  
ma1   -0.341537   0.040571  -8.4182 < 2.2e-16 ***
ma2    0.204242   0.039030   5.2329 1.668e-07 ***
ma3   -0.692204   0.028275 -24.4809 < 2.2e-16 ***
drift  0.187589   0.039109   4.7965 1.614e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC: 2597.526
Mostrar Código
holt <- holt(y = train[,'gsml'], h = 12, initial = 'optimal')
holt$model
Holt's method 

Call:
 holt(y = train[, "gsml"], h = 12, initial = "optimal") 

  Smoothing parameters:
    alpha = 0.9999 
    beta  = 0.6589 

  Initial states:
    l = -49.0075 
    b = -2.4967 

  sigma:  2.4035

     AIC     AICc      BIC 
5429.938 5430.030 5452.384 

Modelos

Mostrar Código
extract_eq(modelo, wrap = T, use_coefs = T)

\[ \begin{alignat}{2} &\text{let}\quad &&y_{t} = \operatorname{y}_{\operatorname{0}} +0.21\operatorname{t} -0.17\operatorname{co2}_{\operatorname{t}} +\eta_{t} \\ &\text{where}\quad &&(1 -0.8\operatorname{B} +0.1\operatorname{B}^{\operatorname{2}} )\ (1 - \operatorname{B}) \eta_{t} \\ & &&= (1 -0.34\operatorname{B} +0.21\operatorname{B}^{\operatorname{2}} -0.69\operatorname{B}^{\operatorname{3}} )\ \varepsilon_{t} \\ &\text{where}\quad &&\varepsilon_{t} \sim{WN(0, \sigma^{2})} \end{alignat} \]

\[( 1 \ - \ 0.79B \ + \ 0.09B^2) \ \cdot \ (1 - B) \ \cdot \ (y_t - 0.19t) \\ = (1 \ - \ 0.34B \ + \ 0.2B^2 \ - \ 0.69B^3) \ \cdot \ \epsilon_t\]

Supuestos

Mostrar Código
mod %>%
  checkresiduals(theme = theme_bw(), test = F)

Mostrar Código
mod %>%
  residuals() %>%
  shapiro.test()

    Shapiro-Wilk normality test

data:  .
W = 0.99855, p-value = 0.8733
Mostrar Código
mod %>%
  residuals() %>%
  Box.test(type = 'Ljung-Box')

    Box-Ljung test

data:  .
X-squared = 0.028529, df = 1, p-value = 0.8659

Outliers

Mostrar Código
tso(y = train[,'gsml'], xreg = train[,'co2'], tsmethod = 'auto.arima', args.tsmethod = list(stepwise = F, approximation = F))
Series:  
Regression with ARIMA(2,1,3) errors 

Coefficients:
         ar1      ar2      ma1     ma2      ma3   drift     xreg
      0.7971  -0.0955  -0.3391  0.2081  -0.6935  0.2085  -0.1712
s.e.  0.0550   0.0549   0.0402  0.0387   0.0282  0.0415   0.0877

sigma^2 = 2.993:  log likelihood = -1289.86
AIC=2595.73   AICc=2595.95   BIC=2631.63

No outliers were detected.

Predicciones

Mostrar Código
prediccion <- forecast(object = modelo, h = 12, level = 95, xreg = test[,'co2'])

serie <- prediccion %>%
  as_tibble() %>%
  ts(start = c(2013, 1), end = c(2013, 12), frequency = 12) %>%
  {cbind(data[, 'gsml'], .)}

colnames(serie) <- c('gsml', 'predicho', 'lwr', 'upr')

serie %>%
  dygraph(height = 350) %>%
  dySeries('gsml', label = 'Observada') %>% 
  dySeries(c('lwr', 'predicho', 'upr'), label = 'Predicho') %>%
  dyRangeSelector()
Mostrar Código
accuracy(prediccion, test[,'gsml'])[,c(2,5)]
                 RMSE     MAPE
Training set 1.719463      Inf
Test set     9.180928 12.04332
Mostrar Código
serie <- holt %>%
  as_tibble() %>%
  select(1, 4, 5) %>%
  rename(predico = 1, lwr = 2, upr = 3) %>%
  ts(start = c(2013, 1), end = c(2013, 12), frequency = 12) %>%
  {cbind(data[, 'gsml'], .)}

colnames(serie) <- c('gsml', 'predicho', 'lwr', 'upr')

serie %>%
  dygraph(height = 350) %>%
  dySeries('gsml', label = 'Observada') %>% 
  dySeries(c('lwr', 'predicho', 'upr'), label = 'Predicho') %>%
  dyRangeSelector()
Mostrar Código
accuracy(holt, test[,'gsml'])[,c(2,5)]
                 RMSE     MAPE
Training set 2.396217      Inf
Test set     4.513079 5.701727

Gracias